Fin most abundant taxa across sites:

sapply(c(biobakery_function, alpha, beta, taxa, norm, heat, varia), source, chdir = TRUE)

Data loading:

input <- list(metpahlan_ps = "~/Documents/GitHub/oral-microbiome/data/processed/metaphlan_ps.RDS",
              motus_ps = "~/Documents/GitHub/oral-microbiome/data/processed/motus/physeq.RDS",
              pathway_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_pw_ps.RDS",
              kegg_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_kegg_ps.RDS",
              l4c_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_l4c_ps.RDS",
              go_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_go_ps.RDS",
              infogo_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_infogo_ps.RDS",
              metac_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_metac_ps.RDS",
              pfam_ps = "~/Documents/GitHub/oral-microbiome/data/processed/humann_pfam_ps.RDS")

figures = "~/Documents/GitHub/oral-microbiome/figures"
dir.create(figures)
## Warning in dir.create(figures): '/Users/fconstan/Documents/GitHub/oral-
## microbiome/figures' already exists
purrr::map(input, 
           readRDS) -> objects
sample_data(objects$metpahlan_ps$physeq)$Site_Health <- paste0(sample_data(objects$metpahlan_ps$physeq)$Oral_Site,
                                                               "_",
                                                               sample_data(objects$metpahlan_ps$physeq)$Health_status)

Fig1 Identification of predominant streptococcus species at all sites

This figure is good, but with the approach we decided focusing solely on streptococcus species We should change it a bit so that it shows top 10 streptococcus species across all three sites

Identification of predominant streptococcus species at all sites

Fig1A: predominant streptococcus species across all samples in health

objects$metpahlan_ps$physeq %>%
  subset_samples(Health_status == "Healthy control") -> tmp
# subset_samples(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) -> tmp

tmp %>%
  subset_taxa(Genus %in% c("Streptococcus")) %>%
  physeq_most_abundant(physeq = .,
                       group_var = "Oral_Site",
                       ntax = 10,
                       tax_level = "Species") -> ma_saliva_tongue_strep
tmp %>%
  transform_sample_counts(function(x) x/sum(x) * 100) %>%
  phyloseq::subset_taxa(Species %in% ma_saliva_tongue_strep) %>%
  phyloseq_ampvis_heatmap(transform = FALSE,
                          group_by = "Subject",#treatment
                          facet_by = c("Oral_Site", "Health_status", "Subject"),
                          ntax = 100,
                          tax_aggregate = "Species",
                          tax_add = NULL) -> p
## Warning: There are only 16 taxa, showing all
p + facet_grid(. ~ Oral_Site + Health_status, scales = "free", space = "free") + 
  scale_fill_viridis_c(breaks = c(0,  0.01, 1, 10, 25, 50, 75, 100), 
                       labels = c(0,  0.01, 1, 10, 25,50, 75, 100), 
                       trans = scales::pseudo_log_trans(sigma = 0.1))  + 
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6)) + 
  theme(axis.text.y = element_text(angle = 0,  size = 8))-> Fig1

Fig1

Fig1 %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))


Fig1 %>%
  ggsave(path = figures,
         plot = .,
         filename = "Fig1.png",
         units = "cm", width = 20, height = 10)
objects$metpahlan_ps$physeq %>%
  subset_samples(Health_status == "Healthy control") %>% 
  subset_samples(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) -> tmp

tmp %>%
  subset_taxa(Genus %in% c("Streptococcus")) %>%
  physeq_most_abundant(physeq = .,
                       group_var = "Oral_Site",
                       ntax = 10,
                       tax_level = "Species") -> ma_saliva_tongue_strep
ma_saliva_tongue_strep
##  [1] "Streptococcus_australis"     "Streptococcus_cristatus"    
##  [3] "Streptococcus_gordonii"      "Streptococcus_infantis"     
##  [5] "Streptococcus_mitis"         "Streptococcus_mutans"       
##  [7] "Streptococcus_oralis"        "Streptococcus_parasanguinis"
##  [9] "Streptococcus_salivarius"    "Streptococcus_sanguinis"    
## [11] "Streptococcus_sp_A12"        "Streptococcus_sp_F0442"     
## [13] "Streptococcus_sp_HMSC071D03"
ma_saliva_tongue_strep_sel <- c("Streptococcus_parasanguinis",
                                "Streptococcus_salivarius",
                                "Streptococcus_infantis")
objects$metpahlan_ps$physeq %>%
  subset_samples(Health_status == "Healthy control") %>% 
  subset_samples(Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) -> tmp

tmp %>%
  subset_taxa(Genus %in% c("Streptococcus")) %>%
  physeq_most_abundant(physeq = .,
                       group_var = "Oral_Site",
                       ntax = 10,
                       tax_level = "Species") -> ma_saliva_plaque_strep
ma_saliva_plaque_strep
##  [1] "Streptococcus_anginosus_group"   "Streptococcus_australis"        
##  [3] "Streptococcus_cristatus"         "Streptococcus_gordonii"         
##  [5] "Streptococcus_infantis"          "Streptococcus_mitis"            
##  [7] "Streptococcus_mutans"            "Streptococcus_oralis"           
##  [9] "Streptococcus_parasanguinis"     "Streptococcus_salivarius"       
## [11] "Streptococcus_sanguinis"         "Streptococcus_sobrinus"         
## [13] "Streptococcus_sp_oral_taxon_056"
ma_saliva_plaque_strep_sel <- c("Streptococcus_oralis",
                                "Streptococcus_sanguinis",
                                "Streptococcus_gordonii")

Fig2 Overall bacterial activity of selected species

Load pathway abundance table and have a look:

objects$pathway_ps$physeq  %>%
  humann2_species_contribution(meta_data_var = c("Subject", "sample_Sample" ,"Type", "Oral_Site", "Health_status")) %>%
  separate(Feature, into = c("code", "name"), sep = ": ", remove = FALSE) %>%
  dplyr::filter(!is.na(name))  -> pwy
## Warning in speedyseq::psmelt(.): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 720 rows [1029,
## 1903, 1961, 2338, 2461, 2854, 3562, 3891, 4578, 5649, 5690, 5849, 5854, 6607,
## 6843, 7028, 7077, 7309, 7819, 8576, ...].
pwy %>%
  distinct(name, .keep_all = TRUE) %>%
  DT::datatable()

Fig2A: box plot species saliva vs. tongue

ma_saliva_tongue_strep_sel
## [1] "Streptococcus_parasanguinis" "Streptococcus_salivarius"   
## [3] "Streptococcus_infantis"
humann2_RNA_DNA_ratio_plot <- function(df,
                                       rm_un_sp = TRUE,
                                       x_plot,
                                       y_plot,
                                       color,
                                       fill,
                                       shape = NULL,
                                       filter_feature = FALSE,
                                       filter_genus = FALSE,
                                       filter_species = FALSE,
                                       only_pos = TRUE, # dplyr::filter(DNA  > 0 , RNA > 0)
                                       facet_formula = FALSE,#". ~ Health_status"
                                       export_legend = FALSE,
                                       box_width = 0.8){ #  
  
  require(tidyverse); require(ggConvexHull)
  
  if(rm_un_sp == TRUE){
    df %>%
      dplyr::filter(Species != "unclassified") -> df
  }
  
  if(filter_feature != FALSE){
    df %>%
      dplyr::filter(Feature %in% filter_feature) -> df
  }
  
  if(filter_genus != FALSE){
    df %>%
      dplyr::filter(Genus %in% filter_genus) -> df
  }
  if(filter_species != FALSE){
    df %>%
      dplyr::filter(Species %in% filter_species) -> df
  }
  
  if(only_pos == TRUE){
    df %>%
      dplyr::filter(DNA  > 0 , RNA > 0) -> df
  }
  
  df %>%
    ggplot(aes_string(x = x_plot, y = y_plot, color = color, fill = fill, shape = shape)) + 
    #geom_boxplot(outlier.colour = NA, alpha = 0.2,
                 # position = position_dodge(width=box_width)) +
    # geom_jitter(size=1, alpha=0.2) +
    # ggbeeswarm::geom_beeswarm(size=1, alpha=0.2) +
    # geom_violin(size=0.05, alpha=0.2,
    #             position = position_dodge(width=box_width)) +
    geom_hline(yintercept = 0,
               col = "red",
               linetype = "dotted",
               size = 0.5) +
    theme_light() +
    guides(fill=guide_legend(ncol=1)) +
    theme(axis.title.x = element_blank()) +
    # ylab(paste0("Relative Abundance (Top ",n," (RNA)) \n")) + #scale_y_continuous(trans='sqrt') + 
    # theme(legend.text = element_text(size= 6)) + 
    # theme(axis.text.x = element_text(size = 4)) +
    theme(legend.key.size = unit(0.2,"cm")) -> p
  
  if(facet_formula != FALSE){
    p + facet_grid(as.formula(facet_formula), drop=TRUE,
                   scale="fixed",space="free_x") + xlab(label = NULL) -> p
  }
  
  if(export_legend == TRUE){
    
    out <- list("legend" = p %>% ggpubr::get_legend() %>% ggpubr::as_ggplot(),
                "plot" = p  + ggpubr::rotate_x_text(90) + 
                  ggpubr::rotate() + theme(legend.position = "none"))
    
  }else{
    out <- p  + ggpubr::rotate_x_text(90) + 
      ggpubr::rotate() + theme(legend.position = "none")
  }
  return(out)
}
pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # filter(RNA_DNA > mean(RNA_DNA, na.rm = TRUE)) %>%
  # dplyr::filter(grepl("ose",name)) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "Species", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_tongue_strep_sel,
                             facet_formula = ". ~ .",
                             export_legend = TRUE) -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
# plot$legend
# 
# plot$plot
ggpubr::get_palette("npg", 3)
## [1] "#E64B35FF" "#4DBBD5FF" "#00A087FF"
plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

# plot$plot + geom_violin(fill = "transparent") + ggforce::geom_sina(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot + geom_violin(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#4DBBD5FF"))  + theme(legend.position = "top") -> fig2a

fig2a 

fig2a %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

fig2a %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig2a.png",
         units = "cm", width = 20, height = 10)
pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # filter(RNA_DNA > mean(RNA_DNA, na.rm = TRUE)) %>%
  # dplyr::filter(grepl("ose",name)) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "Subject", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_tongue_strep_sel,
                             facet_formula = "Oral_Site ~ Species",
                             export_legend = TRUE) -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
# plot$legend
# 
# plot$plot
ggpubr::get_palette("npg", 3)
## [1] "#E64B35FF" "#4DBBD5FF" "#00A087FF"
plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

# plot$plot + geom_violin(fill = "transparent") + ggforce::geom_sina(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot + geom_violin(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#4DBBD5FF"))  + theme(legend.position = "top") -> fig2a

# fig2a 

fig2a %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

# fig2a %>%
#   ggsave(path = figures,
#          plot = .,
#          filename = "fig2a.png",
#          units = "cm", width = 20, height = 10)
plot$plot$data %>% 
  ggpubr::compare_means(RNA_DNA ~ Oral_Site,
                        group.by = c("Species"),
                        data = .,
                        method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05)  %>%
  DT::datatable()

Fig2C: box plot species saliva vs. plaque

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
  # filter(RNA_DNA > mean(RNA_DNA, na.rm = TRUE)) %>%
  # dplyr::filter(grepl("ose",name)) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "Species", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_plaque_strep_sel,
                             facet_formula = ". ~ .",
                             export_legend = TRUE) -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

ggpubr::get_palette("npg", 3)
## [1] "#E64B35FF" "#4DBBD5FF" "#00A087FF"
plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot + geom_violin(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#00A087FF"))  + theme(legend.position = "top") -> fig2b

fig2b

fig2b %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))


fig2b %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig2b.png",
         units = "cm", width = 20, height = 10)
plot$plot$data %>% 
  ggpubr::compare_means(RNA_DNA ~ Oral_Site,
                        group.by = c("Species"),
                        data = .,
                        method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05)  %>%
  DT::datatable()

Fig 3. Pathway activity of selected species

Fig3A: significant pathway expression by candidate species saliva versus tongue in health

redo main plot + stats

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_tongue_strep_sel,
                             facet_formula = " . ~ Health_status + Species",
                             export_legend = TRUE) -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

# plot$plot
plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Oral_Site,
                        group.by = c("Species","Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name",
                             y_plot = "log2(RNA_DNA)",
                             color = "Oral_Site",
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = "Streptococcus_infantis",
                             filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Health_status ",
                             export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot +     geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#4DBBD5FF"))  + theme(legend.position = "top") -> fig3a

fig3a

fig3a %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

fig3a %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig3a.png",
         units = "cm", width = 20, height = 10)
pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = "Streptococcus_parasanguinis",
                             filter_feature = diff_signif %>% filter(Species == "Streptococcus_parasanguinis") %>% pull(Feature),
                             facet_formula = " . ~ Health_status ",
                             export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot + geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#4DBBD5FF")) + theme(legend.position = "top") -> fig3b

fig3b

fig3b %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))


fig3b %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig3b.png",
         units = "cm", width = 20, height = 10)
pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = "Streptococcus_salivarius",
                             filter_feature = diff_signif %>% filter(Species == "Streptococcus_salivarius") %>% pull(Feature),
                             facet_formula = " . ~ Health_status ",
                             export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot +     geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp


ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#4DBBD5FF"))  + theme(legend.position = "top")  -> fig3c

fig3c

fig3c %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

fig3c %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig3c.png",
         units = "cm", width = 20, height = 10)

Fig3B: significant pathway expression by candidate species saliva versus plaque in health

redo main plot + stats

pwy %>%
  dplyr::filter(Health_status == "Healthy control",
                Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = ma_saliva_plaque_strep_sel,
                             facet_formula = " . ~ Health_status + Species",
                             export_legend = TRUE) -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

# plot$plot
plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Oral_Site,
                        group.by = c("Species","Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name",
                             y_plot = "log2(RNA_DNA)",
                             color = "Oral_Site",
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = "Streptococcus_oralis",
                             filter_feature = diff_signif %>% filter(Species == "Streptococcus_oralis") %>% pull(Feature),
                             facet_formula = " . ~ Health_status ",
                             export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot +     geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#00A087FF"))  + theme(legend.position = "top") -> fig4a

fig4a

fig4a %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

fig4a %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig4a.png",
         units = "cm", width = 20, height = 10)
pwy %>%
  dplyr::filter(Health_status == "Healthy control") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Oral_Site", 
                             fill = "Oral_Site",
                             shape = NULL,
                             filter_species = "Streptococcus_sanguinis",
                             filter_feature = diff_signif %>% filter(Species == "Streptococcus_sanguinis") %>% pull(Feature),
                             facet_formula = " . ~ Health_status ",
                             export_legend = TRUE) -> plot


plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot +     geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#E64B35FF", "#00A087FF"))  + theme(legend.position = "top") -> fig4b

fig4b

fig4b %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

fig4b %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig4b.png",
         units = "cm", width = 20, height = 10)
# pwy %>%
#   dplyr::filter(Health_status == "Healthy control") %>%
#      dplyr::filter(Oral_Site %in% c("SALIVA", "SUBGINGIVAL_PLAQUE")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Oral_Site", 
#                        fill = "Oral_Site",
#                        shape = NULL,
#                        filter_species = "Streptococcus_gordonii",
#                        filter_feature = diff_signif %>% filter(Species == "Streptococcus_gordonii") %>% pull(Feature),
#                        facet_formula = " . ~ Health_status ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Fig 5. Overall bacterial activity of selected species in health versus diseases

Fig5A: box plot species in health and disease at each site

strep_sel <- c(ma_saliva_plaque_strep_sel, ma_saliva_tongue_strep_sel)
pwy %>%
  # dplyr::filter(Health_status == "Healthy control",
  # Oral_Site %in% c("SALIVA", "TONGUE_BIOFILM")) %>%
  # filter(RNA_DNA > mean(RNA_DNA, na.rm = TRUE)) %>%
  # dplyr::filter(grepl("ose",name)) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "Species", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = strep_sel,
                             facet_formula = ". ~ Oral_Site",
                             export_legend = TRUE) -> plot
## Warning in if (filter_species != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot + ggpubr::rotate()

plot$plot$data %>% 
  ggpubr::compare_means(RNA_DNA ~ Health_status,
                        group.by = c("Oral_Site", "Species"),
                        data = .,
                        method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05)  %>%
  DT::datatable()
ggpubr::get_palette("jco", 3)
## [1] "#0073C2FF" "#EFC000FF" "#868686FF"
plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot +     geom_violin(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#EFC000FF", "#868686FF"))  + theme(legend.position = "top") -> fig5

fig5

fig5 %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))


fig5 %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig5.png",
         units = "cm", width = 20, height = 10)

Fig 6 Pathway activity of selected species at each site in health versus disease

Fig6A-F: significant pathway expression by candidate species at each site in health and disease

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  # dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             # filter_species = "Streptococcus_infantis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species + Oral_Site ",
                             export_legend = TRUE) -> plot
plot$legend

# plot$plot

A - SALIVA

strep_sel <- c(ma_saliva_plaque_strep_sel, ma_saliva_tongue_strep_sel)
strep_sel
## [1] "Streptococcus_oralis"        "Streptococcus_sanguinis"    
## [3] "Streptococcus_gordonii"      "Streptococcus_parasanguinis"
## [5] "Streptococcus_salivarius"    "Streptococcus_infantis"

Streptococcus_infantis

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_infantis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
# plot$legend
# plot$plot
plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") -> diff

diff %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_infantis") %>%
#      dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Streptococcus_parasanguinis

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_parasanguinis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
# plot$legend
# plot$plot
plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
pwy %>%
  dplyr::filter(Species %in% "Streptococcus_parasanguinis") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             # filter_species = "Streptococcus_infantis",
                             filter_feature = diff_signif %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot + geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#EFC000FF", "#868686FF"))  + theme(legend.position = "top") -> fig6a

fig6a  %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))

fig6a %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig6a.png",
         units = "cm", width = 20, height = 10)

Streptococcus_salivarius

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_salivarius",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
# plot$legend
# plot$plot
plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
pwy %>%
  dplyr::filter(Species %in% "Streptococcus_salivarius") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             # filter_species = "Streptococcus_infantis",
                             filter_feature = diff_signif %>% pull(Feature),
                             facet_formula = "  ~ Species ",
                             export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot + geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp

ggpubr::set_palette(tmp, palette = c("#EFC000FF", "#868686FF"))  + theme(legend.position = "top") -> fig6b


fig6b  %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))


fig6b

fig6b %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig6b.png",
         units = "cm", width = 20, height = 10)

Streptococcus_gordonii

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_gordonii",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_gordonii") %>%
#      dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Streptococcus_sanguinis

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_sanguinis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_sanguinis") %>%
#      dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Streptococcus_oralis

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_oralis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
pwy %>%
  dplyr::filter(Species %in% "Streptococcus_oralis") %>%
  dplyr::filter(Oral_Site %in% c("SALIVA")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             # filter_species = "Streptococcus_infantis",
                             filter_feature = diff_signif %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

B - TONGUE_BIOFILM

strep_sel <- c(ma_saliva_tongue_strep_sel)
strep_sel
## [1] "Streptococcus_parasanguinis" "Streptococcus_salivarius"   
## [3] "Streptococcus_infantis"

Streptococcus_infantis

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_infantis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
plot$legend

plot$plot

plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_infantis") %>%
#      dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
#   humann2_RNA_DNA_ratio_plot(x_plot = "name",
#                        y_plot = "log2(RNA_DNA)",
#                        color = "Health_status",
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif  %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
# 
# 
# plot$legend
# plot$plot

Streptococcus_parasanguinis

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_parasanguinis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
# plot$legend
# plot$plot
plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
pwy %>%
  dplyr::filter(Species %in% "Streptococcus_parasanguinis") %>%
  dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
  humann2_RNA_DNA_ratio_plot(x_plot = "name",
                             y_plot = "log2(RNA_DNA)",
                             color = "Health_status",
                             fill = "Health_status",
                             shape = NULL,
                             # filter_species = "Streptococcus_infantis",
                             filter_feature = diff_signif %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
## Warning in if (filter_feature != FALSE) {: the condition has length > 1 and only
## the first element will be used
plot$legend

plot$plot

plot$plot + geom_point(size  = 0.6, alpha = 0.4, position = position_jitterdodge(dodge.width = 0.8)) + ggpubr::rotate() -> tmp

plot$plot +     geom_boxplot(size=0.05, alpha=0.2) + ggforce::geom_sina(size  = 0.6, alpha = 0.4) + ggpubr::rotate() -> tmp


ggpubr::set_palette(tmp, palette = c("#EFC000FF", "#868686FF"))  + theme(legend.position = "top") -> fig7


fig7  %>% 
  export::graph2ppt(append = TRUE, file = file.path(figures, "figures.pptx"))


fig7 %>%
  ggsave(path = figures,
         plot = .,
         filename = "fig7.png",
         units = "cm", width = 20, height = 10)

Streptococcus_salivarius

pwy %>%
  dplyr::filter(Species %in% "Streptococcus_salivarius") %>%
  dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_salivarius",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
# plot$legend
# plot$plot
plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_salivarius") %>%
#      dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

C - SUBGINGIVAL_PLAQUE

strep_sel <- c(ma_saliva_plaque_strep_sel)
strep_sel
## [1] "Streptococcus_oralis"    "Streptococcus_sanguinis"
## [3] "Streptococcus_gordonii"

Streptococcus_oralis

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SUBGINGIVAL_PLAQUE")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_oralis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
# plot$legend
# plot$plot
plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_infantis") %>%
#      dplyr::filter(Oral_Site %in% c("SUBGINGIVAL_PLAQUE")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>%
#   humann2_RNA_DNA_ratio_plot(x_plot = "name",
#                        y_plot = "log2(RNA_DNA)",
#                        color = "Health_status",
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif  %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
# 
# 
# plot$legend
# plot$plot

Streptococcus_sanguinis

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SUBGINGIVAL_PLAQUE")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_sanguinis",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
# plot$legend
# plot$plot
plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_sanguinis") %>%
#      dplyr::filter(Oral_Site %in% c("SUBGINGIVAL_PLAQUE")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot

Streptococcus_gordonii

pwy %>%
  dplyr::filter(Species %in% strep_sel) %>%
  dplyr::filter(Oral_Site %in% c("SUBGINGIVAL_PLAQUE")) %>%
  # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
  humann2_RNA_DNA_ratio_plot(x_plot = "name", 
                             y_plot = "log2(RNA_DNA)", 
                             color = "Health_status", 
                             fill = "Health_status",
                             shape = NULL,
                             filter_species = "Streptococcus_gordonii",
                             # filter_feature = diff_signif %>% filter(Species == "Streptococcus_infantis") %>% pull(Feature),
                             facet_formula = " . ~ Species ",
                             export_legend = TRUE) -> plot
# plot$legend
# plot$plot
plot$plot$data %>% 
  # mutate(log2RNA_DNA = log2(RNA_DNA)) %>%
  ggpubr::compare_means(RNA_DNA ~  Health_status,
                        group.by = c("Feature"),
                        data = .,
                        # method = "wilcox.test",
                        p.adjust.method = "fdr") %>%
  filter(p.adj < 0.05) -> diff_signif  

diff_signif %>%
  DT::datatable()
# pwy %>%
#    dplyr::filter(Species %in% "Streptococcus_salivarius") %>%
#      dplyr::filter(Oral_Site %in% c("SUBGINGIVAL_PLAQUE")) %>%
#    # dplyr::filter(Oral_Site %in% c("TONGUE_BIOFILM")) %>% 
#   humann2_RNA_DNA_ratio_plot(x_plot = "name", 
#                        y_plot = "log2(RNA_DNA)", 
#                        color = "Health_status", 
#                        fill = "Health_status",
#                        shape = NULL,
#                        # filter_species = "Streptococcus_infantis",
#                        filter_feature = diff_signif %>% pull(Feature),
#                        facet_formula = " . ~ Species ",
#                        export_legend = TRUE) -> plot
#   
# 
# plot$legend
# plot$plot
sessionInfo()
## R version 4.0.5 (2021-03-31)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggConvexHull_0.1.0 gdtools_0.2.3      ampvis2_2.6.4      reshape2_1.4.4    
##  [5] scales_1.1.1       phyloseq_1.34.0    forcats_0.5.1      stringr_1.4.0     
##  [9] dplyr_1.0.7        purrr_0.3.4        readr_1.4.0        tidyr_1.1.3       
## [13] tibble_3.1.2       ggplot2_3.3.4      tidyverse_1.3.1   
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1            uuid_0.1-4              backports_1.2.1        
##   [4] systemfonts_1.0.2       plyr_1.8.6              igraph_1.2.6           
##   [7] lazyeval_0.2.2          splines_4.0.5           fantaxtic_0.1.0        
##  [10] crosstalk_1.1.1         digest_0.6.27           foreach_1.5.1          
##  [13] htmltools_0.5.1.1       fansi_0.5.0             magrittr_2.0.1         
##  [16] cluster_2.1.2           doParallel_1.0.16       openxlsx_4.2.4         
##  [19] Biostrings_2.58.0       modelr_0.1.8            officer_0.3.18         
##  [22] colorspace_2.0-2        rvest_1.0.0             ggrepel_0.9.1          
##  [25] haven_2.4.1             xfun_0.24               crayon_1.4.1           
##  [28] jsonlite_1.7.2          survival_3.2-11         iterators_1.0.13       
##  [31] ape_5.5                 glue_1.4.2              polyclip_1.10-0        
##  [34] rvg_0.2.5               gtable_0.3.0            zlibbioc_1.36.0        
##  [37] XVector_0.30.0          car_3.0-10              Rhdf5lib_1.12.1        
##  [40] BiocGenerics_0.36.1     abind_1.4-5             DBI_1.1.1              
##  [43] rstatix_0.7.0           miniUI_0.1.1.1          Rcpp_1.0.6             
##  [46] viridisLite_0.4.0       xtable_1.8-4            foreign_0.8-81         
##  [49] stats4_4.0.5            DT_0.18                 htmlwidgets_1.5.3      
##  [52] httr_1.4.2              RColorBrewer_1.1-2      ellipsis_0.3.2         
##  [55] pkgconfig_2.0.3         farver_2.1.0            sass_0.4.0             
##  [58] dbplyr_2.1.1            utf8_1.2.1              labeling_0.4.2         
##  [61] tidyselect_1.1.1        rlang_0.4.11            manipulateWidget_0.11.0
##  [64] later_1.2.0             munsell_0.5.0           cellranger_1.1.0       
##  [67] tools_4.0.5             cli_2.5.0               generics_0.1.0         
##  [70] statnet.common_4.5.0    ade4_1.7-17             export_0.2.2.9001      
##  [73] broom_0.7.8             evaluate_0.14           biomformat_1.18.0      
##  [76] fastmap_1.1.0           yaml_2.2.1              knitr_1.33             
##  [79] fs_1.5.0                zip_2.2.0               rgl_0.106.8            
##  [82] nlme_3.1-152            mime_0.11               xml2_1.3.2             
##  [85] compiler_4.0.5          rstudioapi_0.13         curl_4.3.2             
##  [88] plotly_4.9.4.1          ggsignif_0.6.2          reprex_2.0.0           
##  [91] tweenr_1.0.2            bslib_0.2.5.1           stringi_1.6.2          
##  [94] highr_0.9               stargazer_5.2.2         lattice_0.20-44        
##  [97] Matrix_1.3-4            ggsci_2.9               vegan_2.5-7            
## [100] microbiome_1.12.0       permute_0.9-5           multtest_2.46.0        
## [103] vctrs_0.3.8             pillar_1.6.1            lifecycle_1.0.0        
## [106] rhdf5filters_1.2.0      jquerylib_0.1.4         cowplot_1.1.1          
## [109] data.table_1.14.0       flextable_0.6.6         httpuv_1.6.1           
## [112] patchwork_1.1.1         R6_2.5.0                promises_1.2.0.1       
## [115] network_1.17.1          rio_0.5.27              IRanges_2.24.1         
## [118] codetools_0.2-18        ggnet_0.1.0             MASS_7.3-54            
## [121] assertthat_0.2.1        rhdf5_2.34.0            withr_2.4.2            
## [124] S4Vectors_0.28.1        mgcv_1.8-36             parallel_4.0.5         
## [127] hms_1.1.0               grid_4.0.5              speedyseq_0.5.3.9001   
## [130] coda_0.19-4             rmarkdown_2.9           carData_3.0-4          
## [133] Rtsne_0.15              ggpubr_0.4.0            ggforce_0.3.3          
## [136] Biobase_2.50.0          shiny_1.6.0             lubridate_1.7.10       
## [139] base64enc_0.1-3